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ABSTRACT 

Transport by meridional flows has significant consequences for stellar evolution, but is difficult to 
capture in global-scale numerical simulations because of the wide range of timescales involved. Stellar 
evolution models therefore usually adopt parameterizations for such transport based on idealized 
laminar or mean-field models. Unfortunately, recent attempts to model this transport in global 
simulations have produced results that are not consistent with any of these idealized models. In 
an effort to explain the discrepancies between global simulations and idealized models, we here use 
three-dimensional local Cartesian simulations of compressible convection to study the efficiency of 
transport by meridional flows below a convection zone in several parameter regimes of relevance to 
the Sun and solar-type stars. In these local simulations we are able to establish the correct ordering 
of dynamical timescales, although the separation of the timescales remains unrealistic. We find that, 
even though the generation of internal waves by convective overshoot produces a high degree of time 
dependence in the meridional flow field, the mean flow has the qualitative behavior predicted by 
laminar, “balanced” models. In particular, we observe a progressive deepening, or “burrowing”, of 
the mean circulation if the local Eddington-Sweet timescale is shorter than the viscous diffusion 
timescale. Such burrowing is a robust prediction of laminar models in this parameter regime, but 
has never been observed in any previous numerical simulation. We argue that previous simulations 
therefore underestimate the transport by meridional flows. 
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1. INTRODUCTION 

1.1. The Role of Meridional Circulations in Transport and Mixing 

The interiors of stars, in common with all rotating fluids, frequently exhibit meridional circulations, and the transport 
of an gular momen t um and chemical species by these circulatio ns has significant consequences for stellar evolution 
le.g.. iMestell 119531 : iChabover fc Zahi]ll992t iPinsonneaultl 1199^ . Within convective zones, transport by meridional 
circulations is generally less important than transport by turbulent convective motions, but within radiative zones 
transport by meridional circ ulations may be a dominant mechanism. Such transport may play a significant role in the 
spin-down of the Sun’s core (iHoward et all 1967t lSDiegellll972l: lClarklll975ll and the depletion of lithium in solar-type 
stars le.g.. ICharbonneau fc Michaudlll988t FCaraud fc Bodenheimeill20in[ l. Meridiona l flows also transport magnetic 
flux, a fact that forms the bas is of certain theories of t he solar magnetic cycle le.g.. iDikoati fc Charbonneaul[l99^ 
and the solar interior rotation (jGough fc Mclntvr^ 1 19981 1. 

In one-dimensional (ID) models of stellar interiors, tran sport by me ridional circulations is sometimes parameterized 
in terms of a “rotationally-induced mixing” parameter ((ZahnI Il992f) derived from idealized laminar or mean-field 
models. Many such models have been proposed, leading to many different prescriptions for transport and mixing, all 
of which assume that the amplitude and radial extent of the circulation, and hence the degree of mixing, depend only 
on the spherically-averaged angular velocity, U(r). Suc h a description is clear ly limited: in the Sun, for example, the 
spherically-averaged angular velocity is almost uniform ( Couvidat et al.ll200.3f ). yet the latitudinal differential rotation 
drives a meridional circulation, as described in Section [1 21 Moreover, the validity of such parameterizations has never 
been confirmed by any 3D numerical model. 

In this paper we use a fully compressible, 3D numerical model to study the behavior of meridional flows driven by 
differential rotation of the convective envelopes of solar-type stars. We compare the results of our simulations with 
the predictions of both laminar and mean-field models and test the assumptions on which those models are based. In 
the present work we consider only non-magnetic processes; the effects of magnetic fields will be addressed in future 
papers. Ultimately, our aim is to construct a better parameterization for the role of meridional flows in transporting 
angular momentum, chemical species, and magnetic fields in stellar interiors. 

The rest of the paper is structured as follows. In Sections II.2III.3I we briefly review the physical mechanisms that 
give rise to meridional flows and discuss their expected behavior in the context of stellar interiors. In Section [L4] we 
summarize the results of recent global numerical simulations of angular momentum transport in the solar interior. Our 
numerical model is described in Section [2l and parameter constraints are outlined in Section [3l In Section |4] we present 
results from four simulations performed in different parameter regimes. Our findings are summarized in Section [S] and 
discussed in relation to the results of other studies. 
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1.2. Driving of Meridional Circulations 

The presence of meridional circulations in stellar interiors can be explained by two complemen tary arguments. The 
first is a generalization of the classical Vogt-Eddington argument (jVogtl lEddingtMll925h . In a rotating star, 
a balance between centrifugal, Coriolis, pressure, and gravitational forces in the meridional plane generally requires 
that temperature is non-constant within each horizontal surface, implying that the radiative heat flux has nonzero 
divergence. Within a radiative (i.e., non-convecting) zone, advection of entropy by a meridional flow is the only 
mechanism that can balance the divergence of the heat flux. A meridional circulation must therefore be present in 
order to maintain local thermal equilibrium and circumvent the so-called “von Zeipel paradox” (e.g., IMestell [19991 

Eor a given internal rotation profile, we can in principal deduce the meridional flow required to maintain the balances 
just described. Although it is tempting to say that the rotation profile “drives” the meridional flow, is it more correct 
to say that the persistence of the rotation profile, on timescales for which meridional force balance and local thermal 
equilibrium are expected to apply, implies the presence of a meridional flow. The transport of angular momentum by 
this meridional flow will feed back onto the rotation profil e, causing it to evolve over time. The classical example of 
this problem considers a star initially in uniform rotation (jSweetlll95fl[ l: in that case the meridional circulation that 
arises is known as the Eddington-Sweet circulation. Advection of angular momentum by this circulation causes the 
star to develop differential rotation on the Eddington-Sweet timescale, 

= (S) 

where 17 is the initial rotation rate, N is the buoyancy frequency, k is the thermal diffusivity, and R is the stellar radius. 
For a solar-type star, this timescale is typically longer than the star’s main-sequence lifetime, and so transport by the 
Eddington-Sweet circulation is usually neglected in stellar evolution models. However, in a differentially rotating star 
the meridional flows implied by the Vogt-Eddington argument can be much stronger than the classical Eddington- 
Sweet circulation, particularly in regions where the angular velocity gradient is large. In that case transport by the 
circulation can be significant. 

In the classical Eddington-Sweet problem it is assumed that the interior of the star is not subject to any internal 
torques arising from viscous, Maxwell, or small-scale Reynolds stresses, and so meridional circulations are the only 
means of transporting angular momentum. In fact, the consideration of such torques provides the second argument 
for the presence of meridional flows. In a quasi-steady state, any torque arising from viscous. Maxwell, or Reynolds 
stresses must be balanced by a Coriolis torque, because neither pressure nor gravity has a mean azimuthal component, 
and a meridional flow must be present to provide this Co riolis torque. T he process by which an applied torque drives 
a meridional flow is often called “gyroscopic pumping” (iMcIntvrel 1200011 by analogy with the Ekman pumping that 
occurs within Ekman layers. One example is the solar convection zone, in which turbulent Reynolds stresses ma intain a 
state of differential rotation by systematically transporting angular momentum from high to low latitudes le.g.. lMieschl 
1200511 . This turbulent Reynolds torque, prograde in low l atitudes and retrogr ade in high latitudes, also gyroscopically 
pumps a meridional circulation, as originally discussed bv iKippenhahnl (jl96.‘ll l. In general, part of the circulation must 
extend beneath the convection zone and into the radiation zone (un less by chance t he vertically-integrated pumping 
torque within the convection zone happens to be exactly zero — see IMcIntvrel (|2007l . §8.2)). 

The two arguments summarized here make slightly different assumptions about the balances of momentum and 
internal energy, and so one or the other may be preferred under different circumstances. The gyroscopic pumping 
argument is in a sense more general, since it does not assume the presence of stable stratification. In practice, the 
two arguments are often complementary; for example, the presence of meridional flows below the solar convection 
zone can also be inferred from the differential rotation of the solar tachocline, using the Vogt-Eddington argument 
(iSoiegel fc Zahnl[l992t iGough fc Mclntvi^ 1 199811 . 

1.3. Predictions from Laminar and Mean-Field Models 

An important property of gyroscopically pumped meridional circulations is their tendency to “burrow” through 
stably stratified regions; that is, the circulation extends progressively deeper over time. In this way, a meridional 
circulation that is gyroscopically pumped in the convective envelope of a solar-type star can burrow into the radiative 
interior, exchanging angular momentum between the two zones. 

This burrowing process was originally studied in connection with the problem of solar spin-down, that is, the the 
gradual extraction of angular momen tum from the solar interior caused by magnetic braking at the Sun’s surface 
(|Schatzrnaiil 119621 : iHoward et al.lll967ll. Such studies usually assume d idealized, laminar conditions, or else relied on 
laboratory analogies (|Sakuraill970HBenton fc Clarklll973 : lClarldri975ll . These studies neglected any m ean eff e cts ar ising 
from waves, instabilities, or turbulence. The first description of the burrowing process was given bv iClarfl (|1973[1 . for 
linear perturbations to a state of uniform rotation in a cylinder of stably stratified fluid. Assuming a balance of 
meridional forces (i.e., hydrostatic and cyclostrophic balance), as well as local thermal equilibrium, he showed that a 
change in the rotation of the boundaries of the cylinder drove a meridional circulation within a boundary layer whose 
thickness grew “hyperdiffusively” as Advection of angular momentum by the meridional circulation caused the 
angular velocity perturbation imposed at the boundary to propagate into the cylinder at the same rate. 

Burrowing of meridional circulations is now known t o be a robust feature of models in which meridional adv ection 
dominates the transport of angular momentum (e.g., iHavnes et al.l 119911 : ISoiegel fc ZahnI [19921: lElliottI [l997[) . For 
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fluids that are “heavily stratified”, meaning that the buoyancy frequency N far exceeds the rotation rate burrowing 
requires the presence of a “thermal relaxation” mechanism that mitigates the effect of the stratification. For this 
reason, in stellar interiors the rate of burrowing is dependent o n the rate of radiative diffusion, and so the bur r owing 
pro cess in this conte xt is sometimes called “radiative spreading” (|SDieefel fc Zahnl[l992ti . Following iHavnes et al.l (j 199 11 1 
and IMcIntvrel (I2002f) we prefer to call it “burrowing” because of the circulation’s tendency to extend downward, rather 
than upward, in a fluid with a finite density scale height. 

Assuming that the results of laminar models can be applied directly to stellar interiors, the time required for 
meridional circulations to burrow across an entire stellar radiation zone, and thereby communicate the spin-down of 
the surface all the way to the center, is of the order of the Eddington-Sweet time, tES, given by Equation o, where 
now R is the radius of the radiation zone. Thes e spin-down circulations are closely analogous to the Eddington-Sweet 
circulations described in the previous section (|Clarkl[T97^ . Because the timescale Ies is typically longer than the 
main-sequence lifetime of a solar-type star, spin-down by meridional circulations is not expected to operate all t he way 
to th e center of the s tar. I n the solar interior, the Eddington-Sweet timescale is currently ~ 10^^ years (e.g.. iGoughI 
[2001 . and so iClarkI (Il975f) predicted that the solar rotation rate increases significantly with depth in the radiation 
zone. 

This picture of solar spin-down by burrowing meridional circulations was challenged by the advent of helioseis¬ 
mology, which revealed that the solar radiation zone rotates uniformly in both radius and latitude. This uniform 
rotation is incompatib le with any model in which angular momentum transport occurs only by meridional advection 
([Spiegel &: Zahnlll992ll . implying that other transpor t processes are operat i iig in the solar rad iation zone. Different 
authors have s uggested that an i sotropic turbulen ce ([Spiegel fc Zahniri992l : iChabover fc Zahnlll992H. internal grav¬ 
ity waves fe.g.. lSchatzmanlll993HZahn et ai1ll997ll. or primordial magnetic fields (iCharbonneau fc MacGregodll993l : 
Riidiger fc Kitchatinovlll997tlGough fc Mclntvrelll99^ might be responsible. Various parameteri zations for these pro¬ 
cesses have subsequently been incorporated into one-dimensional stellar evolution models (^e.g.. iCharbonnel fc Talon! 
120051: lEggenberger et aLll2005f ). 


1.4. Results from Global Numerical Simulations 

Recently, attempts have been made to model an gular momen t um transport in t he solar interior using global-scale, 
two- and three-dimensional numerical simulations ([Rogersll2nill : l^un et al.ll2011l : IStrugarek et ahll^llh . These sim¬ 
ulations include the effects of convective turbulence, gravity waves, and magnetic fields self-consistently, and can 
therefore be used to test the predictions of the laminar and mean-held models described above. Unfortunately, these 
studies arrive at rather differe nt conclusions, and none of them lends support to any of the pictures of angular mo¬ 
mentum transport just listed. iBrun et al.l (120111) hnd th at angular momentum transport below the convection zone 
is dominated by viscous stresses, whereas iRogerj (120111) hnds that transport by meridional advection and viscous 
stresses cance l, with apparently no net exchange of angular momentum between the convection and radiation zones. 
Rogerj (120111) als o hnd s that the presence of a magnetic held has little effect on angular momentum transport, whereas 
Striigarek et al.l ([201 If ) hnd that a global-scale interior magnetic held dominates the angular momentum transport, 
but tends to produce nonuniform rotation within the radiation zone. 

None of these simulations exhibit the burrowing of meridional hows predicted by laminar models; in fact, such 
burrowing has never been observed in any self-consistent, three-dimensional numerical simulation. However, the pa¬ 
rameter regime in which burrow ing is expected to occur is rather difficult to achieve in numerical simulations. As 
originally noted bv iClarfl ([1973| ). the transport of angular momentum by meridional circulations acts in competition 
with transport by viscosity. Since the timescale for transport by meridional circulation is very long, of order IeSi 
the transport will be do minated by viscous s t resses unle ss the Prandtl number is very smal l (further detail is given 
in Section [3]). Recently. iGaraud fc Brummelll (I2008D and iGaraud fc Acevedo- Arreguinl (I2009D have studied the pene¬ 
tration of meridional hows into stellar radiative zones using laminar, axisymmetric, steady-state models. They hnd 
that the ratio of ^es and the viscous diffusion timescale plays an important role in determining both the magnitude 
and structure of meridional hows. These results have not yet been conhrmed by any self-consistent, three-dimensional 
numerical model. 

In an attempt to understand the discrepancies between the results of the global-scale numerical simulations, as well 
as their departures from the predictions of earlier models, we here study in greater detail the processes that transport 
angular momentum between the convective envelope and radiative interior of solar-type stars. We use a local Cartesian 
numerical model that incorporates the nonlinear effects of convection and gravity waves, and that allows us to study 
parameter regimes that are not attainable in global-scale simulations. 

2. MODEL 

We use the same code used bv iBrummell et ah! (I2002f ) to study the penetration of convection into a radiation zone. 
This is a fully compressible, pseudo-spectral, /-plane code that solves the ideal gas equations within a Cartesian box 
in a rotating frame. We adopt Cartesian coordinates in which x, y, and z correspond to azimuth, colatitude, and 
depth, respectively. The computational domain, illustrated in Figure [1] is periodic in both horizontal directions, x 
and y. We model a localized region at the interface between the convection zone and radiation zone. Using a local 
Cartesian model, rather than a global spherical model, has the advantage that all of the available computational 
power is devoted to studying the interaction between the two zones. However, our results need to be interpreted in 
the context of the global stellar interior. For simplicity we take the rotation axis to be vertical (H = —De^), which is 
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a rea sonable approximat ion at high latitudes, where the burrowing of meridional flows is expected to be most effective 
(e.g.. lHavnes et al]ll991h . Studies in global models (e.g.. lEllio^ll997f) have shown that burrowing also occurs at lower 
latitudes, and that the direction of burrowing is then roughly parallel to the rotation axis. For this reason, local 
and global studies of meridional circula tions tend to produce results that are qualitatively similar, but differ by a 
“geometrical factor” of order unity ('e.g.. iGaraud fc Bodenheimeiil2010n . 



Figure 1. Illustration of the computational domain. T he w hite and black coloring respectively indicates positive and negative values of 
Ux at time t = to in the simulation described in Section 14.21 The thickness of the convective layer, L, is fixed, but convective overshoot 
causes the shear flow to extend into the radiation zone. 


As in the study of iBrummell et ahl (|2002f) . we choose a vertical profile for the thermal conductivity, k{z), and impose 
a uniform vertical heat flux H at the bottom of the domain. We choose k{z) and H such that the radiative temperature 
gradient is subadiabatic in the lower part of the domain, and superadiabatic (i.e., convectively unstable) in the upper 
part of the domain. This naturally leads to the formation of a lower radiation zone and an upper convection zone. In 
all the computations presented here the box size is 2L x2L x 4L, where L is the thickness of the convection zone. The 
radiation zone is chosen to be three times thicker than the convection zone in order to minimize the influence of the 
bottom boundary on the dynamics. The numerical resolution is typically 100 x 100 x 400; we require greater resolution 
in the vertical direction in order to accurately resolve any boundary layers that form at the interface between the 
convection and radiation zones. The top and bottom boundaries of the domain are both impenetrable and stress-free. 
We impose constant temperature Tq at the top of the domain, z = 0, and a constant heat flux at the bottom. Initially, 
the fluid is at rest and in hydrostatic balance with uniform vertical heat flux throughout, and has pressure po and 
density po at the upper boundary z = 0. The ideal gas equations are nondimensionalized using the thickness of the 
convective layer, L, as the lengthscale, and Ljc as the timescale, where c = y^pojpo is the isothermal sound speed at 
the top of the domain. The temperature, T, pressure, p, and density, p, are nondimensionalized using Tq, po, and po, 
respectively. The dimensionless ideal gas equations then take the form 


dt 


M • V ) M — 2Qpez X u = — Vp -I- gpsz + 2pV • D -|- 


dt 


+ u ■ V ] p = —pV • u 
p = pT 


pT ( ^ -h M • V ) ln(pi/^/p) = V • {k{z)VT) -h ^2pD : D 


( 2 ) 

(3) 

(4) 

(5) 


where 7 = 5/3 is the ratio of specific heats, the constants fl and g are the dimensionless rotation rate and gravitational 
acceleration, and D is the deviatoric rate-of-strain tensor. 


Dij — 


1 dui 

2 dxj 


1 duj 

2 dxi 


1 


'V ■ u6i 


( 6 ) 


We have also introduced the dimensionless dynamic viscosity p and thermal conductivity fc, which are related to their 
dimensional counterparts, p* and k* say, by 


pocL 


pocL ’ 


p = 


and 


(7) 
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where Cp is the specific heat capacity. We take /r to be constant throughout the domain, whereas for k we impose a 
vertical profile of the form 


k{z) 


__ 

1 + exp(20(z — 1)) 


_ ^2 _ 

1 + exp(20(l — z)) ’ 


( 8 ) 


so that fc = fci in the upper layer, z < 1, and fc = /c 2 > fci in the lower layer, z > 1, the change occurring across a 
region of dimensionless thickness ~ 0.1. The bottom of the convection zone is therefore fixed at z = 1, but convective 
motions are able to overshoot into the radiation zone. 

We write the three components of the velocity field as u = {ux,Uy,Uz). In our horizontally-periodic Cartesian 
model, differential rotation corresponds to any x-averaged flow in the x direction, and angular momentum is essentially 
equivalent to azimuthal momentum, pux- Because our computational domain is horizontally symmetric, the Reynolds 
stresses in the convective layer are not able to drive any mean differential rotation. In order to mimic the generation of 
differential rotation in the solar convection zone, we add a volume forcing term to the x component of the momentum 
equation (ED, 

F = X{z,t)p{uT{y,z,t)-Ux). (9) 


The effect of the forcing term is to push the flow toward a prescribe d “target” shear flow u t at a rate A. We have 
chosen to consider a situation analogous to the thought-experiment of lSoiegel fc Zahiil (|1992f) . in which the radiation 
zone is initially in uniform rotation despite the presence of differential rotation in the convection zone. We therefore 
take A to be constant initially, A = Aq, and the target velocity ut is taken to be 


uo(l - z) sin(7ry) 

1 -I- exp(20(z — I)) ’ 


( 10 ) 


which tends to zero at depths z > 1. We thereby suppress differential rotation in the radiation zone and drive 
differential rotation in the convection zone. By suppressing differential rotation in the radiation zone we also prevent 
the burrowing of meridional circulations, which are tied to the differential rotation through the balances described in 
Section 11.21 

In all the computations presented here we take Aq = 217, so that the forcing rate matches the rotation rate, and 
uo = 2kl/TT, so that the shearing timescale is comparable to the rotation period. The Rossby number il~^duT/dy 
is therefore of order unity, and so the d ifferential rotation is somewhat stronger than that observed in the Sun and 
most solar-type stars (e.g.. IReinersll20d^ . By forcing a strong differential rotation, we can study the nonlinear effects 
of finite Rossby number, which were neglected in most of the idealized models reviewed in Section 11.31 Once the 
differential rotation reaches a statistically steady state, at t = to say, the forcing parameters are changed to 


and 


1 -I- exp(20(z — 1)) 

Ut = uo(l — 2:) sin(7ry) 


for t > to 


( 11 ) 

( 12 ) 


so that A tends to zero at depths z > 1. This means that the suppressive forcing is switched off below the convection 
zone for t > to, allowing the differential rotation to propagate into the interior, if the dynamics so dictate. 


3. CHOICE OF PARAMETERS 

The burrowing described in Section 11.31 is expected to occur only under specific parameter conditions, which are 
difficult to achieve in a computational model. The relevant parameter regime can be characterized as a hierarchy of 
timescales in the radiation zone: 

acoustic time <C buoyancy time <C rotation time <C Eddington-Sweet time <C viscous time 
fe.g.. lClai^ll973ll . In terms of our dimensionless parameters, these conditions can be expressed as 

nV 


^ i « 

where JV is the dimensionless buoyancy frequency. 


1 


■C 


217 J kij p 


< 


m/p 


= — 
T 


9^ 


7 


1 dr 

9 dz 


(13) 


(14) 


and IT is a characteristic horizontal lengthscale, which is of order unity in our Cartesian model. We note in particular 
that the viscous time (i.e., the timescale for viscous diffusion across the domain) must exceed the Eddington-Sweet 
time. This condition can be expressed as a constraint on the Prandtl number p,/k 2 within the radiation zone: 


p/k2 


'2ny 
^ ’ /a ) 


(15) 
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This constraint is particularly stringe nt because the second inequal i ty in Equation m implies that the right-hand 
side of Equation (IT^ is <C 1. Following iGaraud fc Acevedo- Arreeuinl (|2009ll we introduce the dimensionless parameter 


2fiy k2 


(16) 


in order to express condition (USD more succinctly as ct <C 1. If thi s condition i s not met then viscou s transport 
of angular momentum dominates the transport by meridional flows (iClarld[19731 : iBenton fc ClarkI 11971 . §3.4), and 
we expec t the burrowing of the merid i onal c irculation to be less efficient. For example, in the laminar, steady-state 
model of IGaraud fc Acevedo- Arreguinl (j2009H . cases with cr > 1 have meridional circ ulations that decay exponen tially 
ben eath the convec tion zone, across a boundary layer similar to that described by iBarcilon fc Pedloskvl (|1967r) (see 
also iLinevkinl 1195511 . 

For any solar-type star at the start of the main-sequence, condition m holds throughout the radiation zone, but as 
the star spins down through magnetic braking, this condition may cease to hold in the most strongly stably stratified 
regions. In the solar interior, for example, condition (1151) holds only within the outer part of the radiation zone 
(jGaraud &: Acevedo- Arreguinll200^ and a very small region at the center. 

Gomputational limitations make it difficult to achieve the “low-sigma” regime described by condition m in a 
numerical model. If the rotation rate 12 and buoyancy frequency N take realistic stellar values, then the Prandtl number 
fj./k 2 must be extremely small, which requires very high spatial resolution. For thi s reason, all global simulations 
of the solar interior have been performed in the “high-sigma” regime, cr > 1 fe.g.. lRogersll201ll : iBrun et al.l 120111 : 
IStrugarek et al.ll2011ll . It is therefore important to test whether the behavior of meridional flows in nonlinear numerical 
simulations depends on a in the same manner as is predicted by laminar models. In order to reach a parameter regime 
with cr < 1, we have chosen here to use a local model, and to adopt more modest values for S2 and N, while still 
preserving the ordering of time s cales given by (1131) . 

The code of iBrummell et akl (I2002D requires the specification of six dimensionless input parameters, which we take 
to be: 

• the values fci and k 2 of the thermal conductivity in the convection and radiation zones; 

• the dynamic viscosity /c; 

• the gravitational acceleration g; 

• the rotation rate 12 ; 

• the heat flux H = k 2 dTjdz at the bottom of the domain. 

For any particular choices of H and ki we can choose values for the remaining parameters in order to satisfy the four 
constraints in (USD; different choices for H and fci correspond to different degrees of compressibility and convective 
overshoot Q In practice, we choose values for H and ki such that the pressure scale height is comparable to the height 
of the domain, and such that convective overshoot extends a distance of order 0.1 below the base of the co nvection 
zone (see Figure [ID, as is thought to be typical for solar-type stars ('e.g.. iBrummell et al.ll2002HRem ^ 12001 . 


4. RESULTS 

In the following sections we present results from four simulations performed in different parameter regimes. The 
first simulation. Case 0, differs from the others in that the conductivity of the upper layer, fci, was chosen to make 
that layer adiabatic, and hence marginally stable to convection. The flow in this simulation therefore remains laminar, 
with no convection or internal wave generation, and so we expect the results to follow the predictions of the laminar 
models summarized in Section o Case 1 has the same parameters as Case 0 for the lower layer, but has a smaller 
conductivity fci in the upper layer, so that the upper layer is convectively unstable. Both Case 0 and Case 1 obey the 
ordering of timescales in (fT^ . and so both have cr < 1; we will refer to this as the “low-sigma regime”. 

Case 2 and Case 3 both have a larger viscosity /i and smaller thermal conductivity k throughout the domain, 
relative to Case 1, such that cr > 1; we will refer to this as the “high-sigma regime”. Case 2 has a larger viscosity and 
conductivity than Case 3, but both have the same Prandtl number n/k, equal to 0.24 in the lower layer, z > 1. In 
both cases the heat flux H was chosen so that the stratification profile, and hence the buoyancy frequency N, matches 
that of Cases 0 and I. Case 2 is only weakly convective, because of its relatively high viscosity, whereas Case 3 is more 
strongly convective, because of its low thermal conductivity. 

In each simulation the dimensionless rotation rate and gravitational acceleration were taken to be 12 = 9.6 x 10“^ 
and g = 0.24 respectively. The other parameters, and relevant timescales, are listed in Table |TJ For reference. Tabled] 
also lists ch aracteristic val ues for the same dimensionless timescales in the solar interior, based on the values at O.7i?0 
reported bv lCoiighl (j2007[ l. To allow the most direct comparison with the numerical simulations, the solar timescales 
are calculated using the lengthscale W = 0.7Rq ~ 4.9 x 10^° cm, and quoted in units of L/c, where L = 1.4 x 10® cm 
is one quarter of the pressure scale height and c = \JvlP — 1-8 x 10^cms“^ is the isothermal sound speed. 


IBrummell et all measure the degree of compressibility and overshoot in terms of the quantities Q = Hjki and S = 


{ Qk2 _ 7 

\ H 7-1, 


)/~ convection simulations presented here have 6 = 0.1 and S = 15. 
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Table 1 

Parameters of each simulation, and corresponding solar parameters 


Case 

lO^fcl 

lO^fca 

lO^/i 

IQ^H 

1 

fV 

1 

/ AT Y 

V2t2/ kolp 

p/2 

nip 

<7 

0 

15.1 

24.1 

0.145 

1.4 

11 

52 

l.OxlO'^ 

7.9x10-1 

0.36 

1 

14.5 

24.1 

0.145 

1.4 

11 

52 

l.OxlO'^ 

7.9x10-1 

0.36 

2 

5.12 

8.53 

2.05 

0.51 

11 

52 

2.9x10^^ 

0.56x10-1 

2.27 

3 

1.93 

3.22 

0.774 

0.19 

11 

52 

7.6x10^^ 

1.5xl0i 

2.27 

Sun 





16 

2400 

4.8x101® 

l.lxlOi® 

0.21 


Note. — The las t fiv e columns give approximate values for the dimensionless 
timescales in Equation 03 and the value of cj, calculated using the time-averaged density 
p and buoyancy frequency N below the overshoot region. 


4.1. Case 0: The Low-Sigma Regime, without Convection 

In this case the upper layer, 0 < z < 1 , is n on-convective, and so we expect the dynamics to follow the predictions of 
the laminar models discussed in Section [1.31 Figure [5] shows the steady azimuthal shear and meridional flow at t = to- 
By this time the flow has reached a steady state with a large-scale azimuthal shear and meridional circulation in the 
upper layer. Both the shear and circulation extend slightly into the lower layer, to an extent that depends on the 
rate at which the target velocity ut tends to zero for z > 1, but for z > 1.3 the flow is exponentially weak. We note 
that the azimuthal shear Ux does not exactly match the target shear flow ut given by Equation (11011 . In the steady 
state, the “residual” forcing Aop(ut — Ux) is balanced primarily by the azimuthal component of the Coriolis force, 
2LlpUy, and this balance determines the strength of the meridional flow within the upper layer. This is an example of 
the process referred to as gyroscopic pumping in Section 11.21 Within the upper layer, the downwelling portion of the 
meridional circulation is stronger than the upwelling portion. This asymmetry is a symptom of the Rossby number 
being of order unity, and is therefore less apparent in the lower layer, where the shear is weaker. 


0.004 

0.002 

0.000 


- 0.002 


-0.004 





Figure 2. Shear flow Ux (left panel) and meridional streamlines (right panel) from Case 0, averaged in x at time t = to. The meridional 
streamlines are drawn as contours of a streamfunction tp, which was computed assuming that V ■ pu = 0. Solid and dashed contours 
indicate anti-clockwise and clockwise circulation respectively. Contour levels for Ux and i/i are cubically spaced to show more detail in the 
radiation zone, where the flows are weakest. 


At t = to the forcing is switched off below z = 1, allowing the flow of the upper layer to propagate into the lower, 
stably stratified layer. The sudden imbalance of forces also generates a spectrum of internal waves. Figure [3] shows 
time averages of the azimuthal shear and meridional flow, averaged over one rotation period, at regular intervals after 
t = to- The time averaging filters out most of the internal wave modes, making the long-time evolution more visible. 
From t = to onwards, the meridional circulation of the upper layer begins to burrow into the lower layer, as indicated by 
the increased vertical extent of the main meridional circulation cell in the lower panel of Figure |31 The mean azimuthal 
shear also propagates into the lower layer, at approximately the same rate, suggesting that the transport of angular 
momentum is the result of advection by the meridional flow, as expected for a laminar flow at these parameter values. 
To verify this, we first use the mass conservation equation ([3]) to write the azimuthal component of the momentum 
equation @ as 

d 

-^{pUx) = -2TLpUy - V • {pUx u) pV'^Ux + F. 


(17) 






























Figure 3. Shear flow and meridional streamlines from Case 0 averaged over one rotation period for successive times t > to ^ using the 
same contour levels as in Figure [2] The time averages were taken over the inter vals Ati, At 2 , Ats, and Af 4 indicated in Figure [4| The 
dashed box in the upper panels indicates the control volume V used in Equation m- 


After integrating over a fixed volume V, and employing the divergence theorem, we obtain 



f dV pUy — J 

' dS • M pUx + p. j 

dV J 

1 dS ■ 'V Ux + 1 

dV J 

f dVF, (18) 

V 

Coriolis 

inertial 

viscous forcing 


where dV represents the boundary of V, and where dS is the area element directed outward. In order to verify that 
the burrowing meridional circulation seen in the lower panel of Figure |3] is responsible for the propagation of the shear 
flow into the interior, we have computed each integral in Equation (1181) for the control volume V indicated by the 
dashed box in the upper panel of Figured The result, after taking a time average over one rotation period, is plotted 
in Figure 01 The left panel of Figured shows that, whereas the azimuthal momentum in the upper layer adjusts to a 
new equilibrium after only a few rotation periods, the lower layer adjusts on a much longer timescale. The right panel 
shows that the long-time transport of azimuthal momentum into the lower layer can be attributed almost entirely 
to the azimuthal Coriolis force arising from the mean meridional flow, which corresponds to advection of angular 
momentum viewed in our rotating frame. 

As a further illustration of the relative contributions of the different processes to the net transport of angular 
momentum, in Figure!^ we present vertical cross-sections of each term on the right-hand side of Equation after 
averaging in azimuth and over the time interval Ati indicated in Eigure 01 Within the upper layer, the Coriolis 
and forcing terms are approximately in balance, and the strength of the meridional flow is determined by gyroscopic 
pumping. Within the lower layer, the Coriolis term dominates all others, leading to the evolution of the differential 
rotation seen in Figure [3l Within this layer, the strength of the meridional flow is determined by meridional force 
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Figure 4. Left panel: The solid line shows the total azimuthal momentum f dVpUx in the dashed box indicated in Figure [3] averaged 
over one rotation period. For comparison, the dashed line shows the tot al az imuthal momentum in the region above the dashed box. Right 
panel: The four terms contributing to the right-hand side of Equation (Hg, and their total, each averaged over one rotation period. The 
length of one rotation period 27r/Q is indicated on the plot. The total duration plotted corresponds to about half of the Eddington-Sweet 
time. 



Figure 5. Vertical cross-sections of each term on the right-hand side of Equation (I17II . and their total, averaged in azimuth and over the 
time interval Ati indicated in Figure [4| The contour levels are cubically spaced, as indicated by the colorbar on the left. Throughout 
the control volume V, indicated by the dashed box, the Coriolis term provides the dominant contribution to the total angular momentum 
transport. 


balance and local thermal equilibrium, as anticipated by the Vogt-Eddington argument of Section [l.2[ and is therefore 
roughly proportional to the vertical gradient of the angular velocity, i.e., to duxjdz in our Cartesian geometry. This 
explains why the differential rotation and meridional circulation propagate at the same rate, and also why the strength 
of the meridional flow, and hence the rate of burrowing, decays over time. The results in this case are entirely consistent 
with the laminar models described in Section o 


4.2. Case 1: The Low-Sigma Regime, with Convection 

We now study the effects of turbulent convective overshoot and internal waves on the burrowing of meridional flows. 
In Case 1, unlike Case 0, the upper layer, 0 < z < 1, is convectively unstable, and has an rms Reynolds number 
~ 550. The turbulent motions in this layer continually generate internal waves, which propagate into the lower, stably 
stratified layer. Nevertheless, for times t ^ to the time-averaged shear and meridional circulation are confined to the 
upper layer, as in Case 0. Motions in the lower layer, z > I, are suppressed by the forcing in that region. 

At t = to the forcing is switched off in the lower layer. Figure [S] shows the mean shear and meridional flow, averaged 
over two rotation periods, at regular intervals after t = to, using the same contour levels as Figure [31 As in Case 0, the 
shear and meridional circulation both extend progressively deeper into the lower layer over time, though in a much less 
regular fashion than in Case 0. In any snapshot of the meridional flow, the mean circulation is completely dominated 
by internal wave motions, but after averaging over two rotation periods we observe a similar mean circulation to that 
in Figure 131 

In FigurejTlwe show the evolution of the azimuthal momentum in the volume indicated by the dashed box in FigurelHl 
as well as the terms contributing to the right-hand side of Equation (1181) for this volume. (This figure is equivalent 
to Figure |4l in Section 14.11 1 Despite the presence of internal waves, we find that the mean meridional circulation 
dominates the long-time transport of azimuthal momentum into the radiative layer, as in Case 0 for which the flow 
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Figure 6. Shear flow and meridional streamlines from Case 1, averaged over two rotation periods for times t > to. Time averages were 
taken over the intervals indicated in Figure [7| The contour levels are the same as in Figure [S] 


was laminar. 




Figure 7. Same plots as Figure [4| but for Case 1 and averaged over two rotation periods. The scale in the right-hand panel is larger than 
in Figure[4| The inertial contributions to the momentum transport, which are associated with waves and turbulent fluctuations, are larger 
than in Case 0, but the Coriolis force from the mean meridional flow remains dominant in the long term, leading to the overall increase in 
azimuthal momentum shown by the solid line in the left panel. The dimensionless Eddington—Sweet time in this case is 10^. 
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Although there is qualitative agreement between the results of Case 0 and Case 1, we note that they differ in certain 
respects. In particular, the propagation of the shear into the lower layer does not proceed monotonically, as can be seen 
in the left panel of Figured (cf. Figure|4]), indicating that the burrowing process does not act continuously throughout 
the simulation. Perhaps more surprising, however, is that the overall rate of burrowing is faster in Case 1 than in 
Case 0. These differences between Case 0 and Case 1 can be attributed to differences in the profile of within the 
upper layer in the two cases. In Case 0, Ux is strongest close to the upper boundary, z = 0, where the forcing is 
strongest. In Case I, Ua, is highly time-dependent, but generally strongest close to the interface z = 1 (see Figure I!]). 
Therefore the amplitude of Ux at the top of the radiation zone is generally stronger in Case I, but varies significantly 
in time, producing a burrowing that is more rapid overall, but very irregular. The occasional reversals in the sign of 
the momentum transport in Figure [7] (e.g., for 2000 < t ^ 3000) correspond to the periods when Ux at the interface 
is weakest. The variations in Ux in the upper layer are quasi-periodic, with a timescale of several rotation periods, 
and may possibly be the result of an inertial mode trapped within the convection zone. (An inertia-gravity oscillation 
within the radiation zone, on the other hand, would necessarily have a period of less than half the rotation period, 
because N > 217.) In that case the oscillation relies on the artificial local geometry of our model, and would not be 
present in a global model. Finally, we note that during the phases when burrowing does occur, the transport of angular 
momentum throughout the radiation zone is qualitatively similar to that in Case 0. This is illustrated in Figure El in 
which we plot vertical cross-sections of each term on the right-hand side of Equation (II3, averaged in azimuth and 
over the interval Ati indicated in Figure [71 (This figure is equivalent to Figure O in Section |4Tl) Although there is a 
non-zero contribution from the inertial and viscous terms, the main contribution comes from the Coriolis term, and 
the overall transport of angular momentum in the radiation zone resembles that in Case 0. 


0.0001 


0.0000 

- 0.0001 







Figure 8. Same plots as Figure[5l but for Case 1 and averaged over the time interval Ati indicated in Figure[7] As in Case 0, the Coriolis 
term provides the dominant contribution to the total angular momentum transport in the radiation zone. 


4.3. Case 2: The High-Sigma Regime; High Viscosity 

In Case 2 the viscous diffusion timescale is shorter than the Eddington-Sweet timescale, primarily because the 
viscosity is larger than in Cases 0 and 1, and so cr > 1. The rotation rate 17 and stratification profile in the radiation 
zone are the same as in Cases 0 and 1. 

Because of the increased viscosity, in this case the upper layer is only weakly convective, and consequently the 
mean flows are more clearly defined. Most of the kinetic energy in the upper, convective layer is associated with the 
mean azimuthal shear flow, rather than with convective motions. The mean azimuthal shear and meridional flow at 
successive times t > to are plotted in Figure El using the same contour levels as in Figures El and [6l We find that the 
convection zone’s azimuthal shear propagates progressively deeper into the radiation zone, as in Cases 0 and I, but that 
meridional flows within the radiation zone remain confined to a thin layer below the bottom of the convection zone. 
The convection zone’s meridional circulation does not burrow significantly, and instead a weak counter-circulating 
meridional cell is established at the top of the radiation zone (cf. Figures El and [H). 

Figure ITOl shows the equivalent of Figures IH and [71 for Case 2. In this case, we find that the propagation of shear 
into the radiation zone is caused primarily by viscous diffusion. After about one viscous diffusion time, a roughly 
steady state is achieved in which the viscous force is balanced by the Coriolis force. We note that the Coriolis force in 
this case has the opposite sign than in Cases 0 and 1, so transport of azimuthal momentum by the mean meridional 
flow acts to oppose (but not prevent) the propagation of the convection zone’s shear into the radiation zone. This 
is a consequence of the counter-circulating meridional cell visible in Figure El as demonstrated by Figure 1111 which 
shows vertical cross-sections equivalent to those of Figures El and El for Case 2. The lack of burrowing of the meridional 
circulation in this case is consistent with the predictions of laminar models with tr > 1, as discussed in Section El 

4.4. Case 3: The High-Sigma Regime; Low Thermal Conductivity 
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Figure 9. Shear and merid ion al flow from Case 2, averaged over one rotation period for times t > to. Time averages were taken over the 
intervals indicated in Figure fTol Contour levels are the same as in Figures [3] and 0 




Figure 10. Same plots as Figures [4| and 0 but for Case 2. The plots have been averaged over one rotation period. In this case the viscous 
term is dominant, and the Coriolis term takes the opposite sign. The total duration plotted corresponds to about one viscous diffusion 
time. 


Case 3 has the same values of fl, N, and a as Case 2, but the viscosity and thermal conductivity are both smaller 
than in Case 2. As a result the upper layer is more strongly convective, with an rms Reynolds number ~ 60. Although 
this Reynolds number is significantly lower than that of Case 1, we find that internal waves are generated with a larger 
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Figure 11. Same plots as Figures [5] and [S] but for Case 2, and averaged over the time interval Ati indicated in Figure [TOl In this case 
the viscous term causes the convection zone’s azimuthal shear to propagate into the radiation zone. The Coriolis term has the opposite 
sign to the viscous term, and a smaller amplitude overall. 


amplitude in this case than in Case 1. We attribute this to the lower value of the thermal conductivity in this case, 
which reduces the damping of internal waves. In fact, we find that a standing gravity mode is excited in the lower 
layer, which slowly grows in amplitude during the simulation. The existence of this mode is dependent on the artificial 
lower boundary of the computational domain, and so we end the simulation at the point where the amplitude of this 
mode becomes large enough to significantly affect the dynamics (at around t = to + 4000). 

Figure [T^ shows the mean azimuthal shear and meridional flow at successive times t > to, using the same contour 
values as Figures [3l [6l and [9] Despite the presence of internal waves, the mean flow exhibits similar behavior to Case 
2. The meridional circulation driven in the convection zone turns over at only a short depth within the radiation 
zone, and a counter-rotating meridional cell forms beneath, though with a more complicated structure than in Case 
2. Meanwhile, the convection zone’s shear propagates monotonically into the interior. 

Figure [T^ shows the equivalent of Figures |H 0 and [TUI for Case 3. As in Case 2, we find that the propagation of 
shear into the radiation zone is caused primarily by viscous diffusion, although the relative contributions from inertial 
and Coriolis forces are larger than in Case 2. In contrast to Case 2, however, the Coriolis term has the same sign as 
the viscous term, at least for times t < 2500. This is because the convection zone’s meridional circulation, shown in 
Figure [121 manages to burrow a short distance into the radiation zone, as illustrated by Figure [CT which shows the 
equivalent of Figures El El and [TT] for Case 3, averaged ove r the time interval At 2 indicated in Figure [131 A similar 
result was obtained by Garaud fc Acevedo- ArreguiiJ (1200911 in a steady-state model of the solar interior; they found 
that, when cr > 1, meridional circulations that are gyroscopically pumped within the convection zone extend a distance 
of order W/a into the radiation zone, where W is the horizontal lengthscale of the circulation, which in our model is 
of order unity. 


5. SUMMARY AND CONCLUSIONS 

This paper presents the first 3D, self-consistent, and nonlinear study of meridional flows in the parameter regime 
described bv iClarlj (|1973f) . which is the relevant parameter regime for the radiative zones of many solar-type stars. 
In this regime the Eddington-Sweet time is shorter than the viscous time; the ratio of these timescales is determined 
by the dimensionless parameter a defined in Equation (fT6)l . We have considered four separate cases: two in the 
“low-sigma” regime and two in the “high-sigma” regime. Our results indicate that the underlying long-time dynamical 
picture of angular momentum transport predicted by laminar models applies even under more realistic conditions, 
including the presence of overshooting from the neighboring turbulent convection zone, and the internal waves that 
this generates. 

In Cases 0 and 1, which have cr < 1, angular momentum transport is dominated in the long term by advection 
by meridional flows, and the meridional circulation driven in the convection zone burrows (i.e., extends progressively 
downward) into the radiation zone on the Eddington-Sweet timescale, carrying the differential rotation of the convec¬ 
tion zone with it. The burrowing is more irregular in Case 1 than in Case 0, as a result of angular momentum transport 
by turbulence and internal waves. In Cases 2 and 3, which have cr > 1, viscous stresses dominate the transport of 
angular momentum in the radiation zone, although there is also a significant contribution from internal waves in Case 
3. In these two cases the meridional circulation driven in the convection zone extends only a short distance into the 
radiation zone. The differential rotation of the convection zone still propagates into the radiation zone, but by viscous 
diffusion rather than by meridional advection. 

It should be borne in mind that the simulations presented here are at most weakly turbulent, in comparison with real 
stellar convection. It may be that under the more turbulent conditions characteristic of real stellar interiors angular 
momentum transport is dominated by shear-driven turbulence or internal wave breaking, as argued for instance by 
iZahnI (|1992f) . If that is the case then our results may not be applicable to real stars. However, ou r results shoul d 
certainly be applicable to previous global-scale simulations of the solar interior, including those of iRoeeri (1201 Ifl . 
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Figure 12. Shear and meridio nal flow from Case 3, averaged over one rotation period for times t > to. Time averages were taken over 
the intervals indicated in Figure fT^ Contour levels are the same as in Figures [3] 0 and [3 




Figure 13. Same plots as Figures [4| [T] and 1 101 but for Case 3. The plots have been averaged over one rotation period. The viscous term 
is dominant, as in Case 2, until t to + 3000. After this time the inertial contribution from a standing gravity mode becomes comparable 
to the viscous term. 


iBrun et all l|201lD . and iStruearek et ^ (|2011l) . We note that these simulations were all performed in the “high- 
sigma” regime in which we found that angular momentum transport is dominated by visc ous stresses. All of th ese 
global models have cr ~ 200 close to the top of the radiation zone (although in the model of lStruearek et ahl (|2011l l cr 
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Figure 14. Same plots as Figures[5l[8l and llll but for Case 3, and averaged over the time interval At 2 indicated in Figure fT^ As in Case 
2, the viscous term is dominant throughout most of the radiation zone, but in this case the Coriolis and inertial terms also contribute to 
the propagation of the convection zone’s shear within a thin layer at the top of the radiation zone. 


drops to around to 20 deeper within the radiation zone). In each of these simulations it was indeed found that viscous 
stresses contribute at leading order to the transport of angular momentum within the radiation zone. The pattern 
of meridional flows found in these simulations is also rather similar to that which we observe in our simulations with 
O' > 1 (Figures [9] and [12), and no burrowing of the circulation was observed. In this situation we expect the transport 
of an gular momentum betw een the convection and radiation zones to occu r on a viscous timescale. This is indeed 
what iBrun et ^ (|2011h and iStrugarek et al.l (j20lH) observe. iRogersI (l2011h reports that an apparently steady state 
is achieved in which viscous and Coriolis forces balance, and uniform rotation is preserved within the radiation zone. 
Based on our results, we suggest that this “steady state” is actually evolving slowly on a viscous timescale. We are 
now conducting global simulations in the low-sigma regime for comparison, to be presented in a forthcoming paper. 
Using a global model will also allow a more realistic study of the effects of internal waves, avoiding the difficulties 
encountered in Case 3. 

All of the simulations presented here have the same rotation rate and stratification profile, and in each case the 
Prandtl number in the radiation zone is smaller than unity. Yet the strength and depth of the mean meridional 
circulations vary drastically between the four cases. This highlights the danger in modeling stellar interiors with 
numerical simulations that have, for example, realistic rotation and stratification but unrealistic diffusivities. Our 
results suggest that the ordering of dynamical timescales is of greater importance than the exact values of those 
timescales when considering angular momentum transport. Real stellar interiors are characterized by the same ordering 
as the simulations presented here, but with substantially greater separation. 

If our results carry over to realistic stellar parameters then they have significant implications not only for angular 
momentum transport, but also for the transport of chemical elements and magnetic flux by meridional flows. In 
particular, the depth to which chemical elements are carried from the convection zone into the radiation zone will 
depend on the depth to which meridional circulations are able to burrow, and hence on the value of a. 

An important issue not addressed in the present work is the contribution of magnetic fields to angular momentum 
transport in stellar interiors. If magnetic fields act t o suppress differential ro tation in the radiation zone, then the 
burrowing of meridional flows will also be suppressed ([Gough fc Mclntvrelll998[ ) . In that case the role of the magnetic 
field is analogous to that of the forcing in the radiation zone in our model for t < to. Effects of magnetic fields will be 
addressed in future papers. 
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